"""EDA峰值检测算法实现

提供多种SCR（皮肤电反应）峰值检测算法：
- 基础峰值检测：基于scipy.signal.find_peaks的简单检测
- 改进峰值检测：更精细的参数调优和后处理
- NeuroKit2峰值检测：集成NeuroKit2库的检测功能

所有检测函数返回峰值位置和属性信息。
"""

import numpy as np
from scipy.signal import find_peaks, peak_widths, peak_prominences
from scipy.ndimage import gaussian_filter1d
from typing import Tuple, Dict, Any
import logging

logger = logging.getLogger(__name__)

def detect_peaks_basic(eda_filtered: np.ndarray, sampling_rate: float = 4.0) -> Tuple[np.ndarray, Dict[str, Any]]:
    """基础峰值检测算法
    
    使用scipy.signal.find_peaks进行基本的SCR峰值检测。
    
    Args:
        eda_filtered: 预处理后的EDA信号
        sampling_rate: 采样率（Hz）
        
    Returns:
        (peaks, properties): 峰值位置数组和属性字典
    """
    try:
        logger.debug(f"开始基础峰值检测，数据长度: {len(eda_filtered)}")
        
        # 处理NaN值
        valid_mask = ~np.isnan(eda_filtered)
        if np.sum(valid_mask) < 10:
            logger.warning("有效数据点太少，无法进行峰值检测")
            return np.array([]), {'peak_heights': np.array([]), 'widths': np.array([])}
        
        # 基础参数设置
        height_threshold = np.nanpercentile(eda_filtered, 75)  # 75分位数作为高度阈值
        distance_min = max(2, int(sampling_rate * 2))  # 最小间隔2秒
        prominence_threshold = np.nanstd(eda_filtered) * 0.1  # 显著性阈值
        width_min = max(1, int(sampling_rate * 0.5))  # 最小宽度0.5秒
        
        # 峰值检测
        peaks, properties = find_peaks(
            eda_filtered,
            height=height_threshold,
            distance=distance_min,
            prominence=prominence_threshold,
            width=width_min
        )
        
        # 过滤掉NaN位置的峰值
        valid_peaks = []
        for peak in peaks:
            if peak < len(valid_mask) and valid_mask[peak]:
                valid_peaks.append(peak)
        
        peaks = np.array(valid_peaks)
        
        # 重新计算属性
        if len(peaks) > 0:
            peak_heights = eda_filtered[peaks]
            
            # 计算峰值宽度
            try:
                widths_result = peak_widths(eda_filtered, peaks, rel_height=0.5)
                widths = widths_result[0] / sampling_rate  # 转换为秒
            except:
                widths = np.full(len(peaks), 0.5)  # 默认宽度
            
            # 计算显著性
            try:
                prominences = peak_prominences(eda_filtered, peaks)[0]
            except:
                prominences = np.abs(peak_heights - np.nanmean(eda_filtered))
            
            properties = {
                'peak_heights': peak_heights,
                'widths': widths,
                'prominences': prominences,
                'peak_indices': peaks
            }
        else:
            properties = {
                'peak_heights': np.array([]),
                'widths': np.array([]),
                'prominences': np.array([]),
                'peak_indices': np.array([])
            }
        
        logger.debug(f"基础峰值检测完成，检测到 {len(peaks)} 个峰值")
        return peaks, properties
        
    except Exception as e:
        logger.error(f"基础峰值检测失败: {str(e)}")
        return np.array([]), {'peak_heights': np.array([]), 'widths': np.array([])}

def detect_peaks_improved(eda_filtered: np.ndarray, sampling_rate: float = 4.0) -> Tuple[np.ndarray, Dict[str, Any]]:
    """改进的峰值检测算法
    
    在基础检测基础上增加：
    1. 自适应阈值设置
    2. 多尺度检测
    3. 峰值质量评估
    4. 假阳性过滤
    
    Args:
        eda_filtered: 预处理后的EDA信号
        sampling_rate: 采样率（Hz）
        
    Returns:
        (peaks, properties): 峰值位置数组和属性字典
    """
    try:
        logger.debug(f"开始改进峰值检测，数据长度: {len(eda_filtered)}")
        
        # 处理NaN值
        valid_mask = ~np.isnan(eda_filtered)
        if np.sum(valid_mask) < 20:
            logger.warning("有效数据点太少，回退到基础检测")
            return detect_peaks_basic(eda_filtered, sampling_rate)
        
        # 1. 自适应参数设置
        signal_std = np.nanstd(eda_filtered)
        signal_mean = np.nanmean(eda_filtered)
        
        # 动态高度阈值（基于信号特征）
        if signal_std > 0.5:
            height_percentile = 80  # 高变异性信号用更高阈值
        else:
            height_percentile = 70  # 低变异性信号用较低阈值
        
        height_threshold = np.nanpercentile(eda_filtered, height_percentile)
        
        # 自适应距离阈值（基于典型SCR持续时间）
        min_scr_interval = max(int(sampling_rate * 1.5), 3)  # 1.5-3秒
        max_scr_interval = int(sampling_rate * 8)  # 最大8秒
        
        # 自适应显著性阈值
        prominence_threshold = max(signal_std * 0.2, 0.01)
        
        # 自适应宽度阈值
        min_width = max(int(sampling_rate * 0.3), 1)  # 最小0.3秒
        max_width = int(sampling_rate * 5)  # 最大5秒
        
        # 2. 多尺度检测
        # 先用宽松参数检测候选峰值
        candidate_peaks, candidate_props = find_peaks(
            eda_filtered,
            height=np.nanpercentile(eda_filtered, 60),  # 更宽松的高度阈值
            distance=min_scr_interval,
            prominence=prominence_threshold * 0.5,  # 更宽松的显著性
            width=min_width
        )
        
        # 3. 峰值质量评估和筛选
        if len(candidate_peaks) == 0:
            logger.debug("未检测到候选峰值")
            return np.array([]), {'peak_heights': np.array([]), 'widths': np.array([])}
        
        # 计算每个峰值的质量分数
        quality_scores = []
        for peak in candidate_peaks:
            if peak >= len(eda_filtered) or not valid_mask[peak]:
                quality_scores.append(0)
                continue
                
            # 高度分数（相对于局部基线）
            local_window = max(int(sampling_rate * 3), 5)
            start_idx = max(0, peak - local_window)
            end_idx = min(len(eda_filtered), peak + local_window)
            local_baseline = np.nanmin(eda_filtered[start_idx:end_idx])
            height_score = (eda_filtered[peak] - local_baseline) / signal_std if signal_std > 0 else 0
            
            # 形状分数（峰值的尖锐程度）
            if peak > 0 and peak < len(eda_filtered) - 1:
                left_slope = eda_filtered[peak] - eda_filtered[peak-1]
                right_slope = eda_filtered[peak-1] - eda_filtered[peak]
                shape_score = (left_slope + right_slope) / 2 if signal_std > 0 else 0
            else:
                shape_score = 0
            
            # 局部显著性分数
            try:
                prominence = peak_prominences(eda_filtered, [peak])[0][0]
                prominence_score = prominence / signal_std if signal_std > 0 else 0
            except:
                prominence_score = 0
            
            # 综合质量分数
            total_score = height_score * 0.4 + shape_score * 0.3 + prominence_score * 0.3
            quality_scores.append(max(0, total_score))
        
        # 4. 基于质量分数筛选峰值
        quality_scores = np.array(quality_scores)
        if len(quality_scores) > 0:
            quality_threshold = np.percentile(quality_scores, 50)  # 保留质量前50%的峰值
            quality_mask = quality_scores >= quality_threshold
            
            # 同时满足基本条件
            height_mask = eda_filtered[candidate_peaks] >= height_threshold
            final_mask = quality_mask & height_mask
            
            peaks = candidate_peaks[final_mask]
        else:
            peaks = np.array([])
        
        # 5. 计算最终峰值属性
        if len(peaks) > 0:
            peak_heights = eda_filtered[peaks]
            
            # 计算峰值宽度
            try:
                widths_result = peak_widths(eda_filtered, peaks, rel_height=0.5)
                widths = widths_result[0] / sampling_rate
                width_heights = widths_result[1]
                left_ips = widths_result[2]
                right_ips = widths_result[3]
            except:
                widths = np.full(len(peaks), 1.0)
                width_heights = peak_heights * 0.5
                left_ips = peaks - 1
                right_ips = peaks + 1
            
            # 计算显著性
            try:
                prominences = peak_prominences(eda_filtered, peaks)[0]
            except:
                prominences = np.abs(peak_heights - np.nanmean(eda_filtered))
            
            # 计算峰值间隔
            if len(peaks) > 1:
                intervals = np.diff(peaks) / sampling_rate
            else:
                intervals = np.array([])
            
            # 计算上升时间和下降时间
            rise_times = []
            fall_times = []
            
            for i, peak in enumerate(peaks):
                # 上升时间（从局部最小值到峰值）
                search_start = max(0, peak - int(sampling_rate * 3))
                pre_peak_segment = eda_filtered[search_start:peak+1]
                if len(pre_peak_segment) > 1:
                    min_idx = np.nanargmin(pre_peak_segment)
                    rise_time = (peak - (search_start + min_idx)) / sampling_rate
                    rise_times.append(rise_time)
                else:
                    rise_times.append(0)
                
                # 下降时间（从峰值到局部最小值）
                search_end = min(len(eda_filtered), peak + int(sampling_rate * 5))
                post_peak_segment = eda_filtered[peak:search_end]
                if len(post_peak_segment) > 1:
                    min_idx = np.nanargmin(post_peak_segment)
                    fall_time = min_idx / sampling_rate
                    fall_times.append(fall_time)
                else:
                    fall_times.append(0)
            
            properties = {
                'peak_heights': peak_heights,
                'widths': widths,
                'prominences': prominences,
                'peak_indices': peaks,
                'quality_scores': quality_scores[final_mask] if len(quality_scores) > 0 else np.array([]),
                'intervals': intervals,
                'rise_times': np.array(rise_times),
                'fall_times': np.array(fall_times),
                'width_heights': width_heights,
                'left_ips': left_ips,
                'right_ips': right_ips
            }
        else:
            properties = {
                'peak_heights': np.array([]),
                'widths': np.array([]),
                'prominences': np.array([]),
                'peak_indices': np.array([]),
                'quality_scores': np.array([]),
                'intervals': np.array([]),
                'rise_times': np.array([]),
                'fall_times': np.array([])
            }
        
        logger.debug(f"改进峰值检测完成，从 {len(candidate_peaks)} 个候选中筛选出 {len(peaks)} 个峰值")
        return peaks, properties
        
    except Exception as e:
        logger.error(f"改进峰值检测失败: {str(e)}，回退到基础检测")
        return detect_peaks_basic(eda_filtered, sampling_rate)

def detect_peaks_neurokit2(eda_filtered: np.ndarray, sampling_rate: float = 4.0) -> Tuple[np.ndarray, Dict[str, Any]]:
    """NeuroKit2峰值检测插件
    
    使用NeuroKit2库进行SCR峰值检测，如果库不可用则回退到改进算法。
    
    Args:
        eda_filtered: 预处理后的EDA信号
        sampling_rate: 采样率（Hz）
        
    Returns:
        (peaks, properties): 峰值位置数组和属性字典
    """
    try:
        import neurokit2 as nk
        logger.debug("使用NeuroKit2进行峰值检测")
        
        # 使用NeuroKit2检测峰值
        signals, info = nk.eda_peaks(eda_filtered, sampling_rate=sampling_rate)
        
        # 提取峰值位置
        peaks = np.where(signals['SCR_Peaks'] == 1)[0]
        
        if len(peaks) > 0:
            # 获取峰值高度
            peak_heights = eda_filtered[peaks]
            
            # 计算峰值宽度（使用NeuroKit2的结果或自己计算）
            if 'SCR_Width' in signals.columns:
                widths = signals['SCR_Width'].iloc[peaks].values
                # 转换为秒（如果需要）
                widths = widths / sampling_rate if np.max(widths) > 10 else widths
            else:
                # 自己计算宽度
                try:
                    widths_result = peak_widths(eda_filtered, peaks, rel_height=0.5)
                    widths = widths_result[0] / sampling_rate
                except:
                    widths = np.full(len(peaks), 0.5 * sampling_rate)
            
            # 获取其他属性
            if 'SCR_Amplitude' in signals.columns:
                amplitudes = signals['SCR_Amplitude'].iloc[peaks].values
            else:
                amplitudes = peak_heights
            
            # 计算显著性
            try:
                prominences = peak_prominences(eda_filtered, peaks)[0]
            except:
                prominences = amplitudes
            
            properties = {
                'peak_heights': peak_heights,
                'widths': widths,
                'prominences': prominences,
                'amplitudes': amplitudes,
                'peak_indices': peaks,
                'neurokit2_signals': signals,
                'neurokit2_info': info
            }
        else:
            properties = {
                'peak_heights': np.array([]),
                'widths': np.array([]),
                'prominences': np.array([]),
                'amplitudes': np.array([]),
                'peak_indices': np.array([]),
                'neurokit2_signals': signals,
                'neurokit2_info': info
            }
        
        logger.debug(f"NeuroKit2峰值检测完成，检测到 {len(peaks)} 个峰值")
        return peaks, properties
        
    except ImportError:
        logger.info("NeuroKit2未安装，使用改进算法")
        return detect_peaks_improved(eda_filtered, sampling_rate)
    except Exception as e:
        logger.warning(f"NeuroKit2峰值检测失败: {str(e)}，使用改进算法")
        return detect_peaks_improved(eda_filtered, sampling_rate)

def _smooth_for_detection(signal: np.ndarray, sampling_rate: float) -> np.ndarray:
    """为峰值检测优化的信号平滑
    
    Args:
        signal: 输入信号
        sampling_rate: 采样率
        
    Returns:
        平滑后的信号
    """
    try:
        # 使用高斯滤波进行轻微平滑
        sigma = sampling_rate * 0.1  # 0.1秒的标准差
        smoothed = gaussian_filter1d(signal, sigma=sigma)
        return smoothed
    except:
        return signal

def _validate_peaks(peaks: np.ndarray, signal: np.ndarray, sampling_rate: float) -> np.ndarray:
    """验证和过滤峰值
    
    Args:
        peaks: 峰值位置
        signal: 信号数据
        sampling_rate: 采样率
        
    Returns:
        验证后的峰值位置
    """
    if len(peaks) == 0:
        return peaks
    
    valid_peaks = []
    
    for peak in peaks:
        # 检查边界
        if peak <= 0 or peak >= len(signal) - 1:
            continue
            
        # 检查是否为局部最大值
        if (signal[peak] > signal[peak-1] and 
            signal[peak] > signal[peak+1]):
            valid_peaks.append(peak)
    
    return np.array(valid_peaks)